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Abstract 

In this paper a new method for inverting the Laplace transform 
from the real axis is formulated. This method is based on a quadrature 
formula. We assume that the unknown function f(t) is continuous 
with (known) compact support. An adaptive iterative method and an 
adaptive stopping rule, which yield the convergence of the approximate 
solution to f{t), are proposed in this paper. 
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1 Introduction 

Consider the Laplace transform : 

roc 

£/(p) == / e-*f(t)dt = Fip), Rep > 0, (1) 
J o 

where £ : X , b -> L 2 [0, oo , 

Xqm ■= {/ G £ 2 [0,oo) | supp/ C [0,b)}, b>0. (2) 



We assume in ([2]) that / has compact support. This is not a restriction 
practically. Indeed, if lim^oo f(t) = 0, then \ f(t)\ < 5 for t > t$, where 5 > 
is an arbitrary small number. Therefore, one may assume that supp/ C 
[0, , and treat the values of / for t > tg as noise. One may also note that 
if / e ^(0,00), then 

poo pb poo 

F(p):= / f(t)e-P t dt= / f(t)e-P t dt+ / f(t)e^ t dt:=F 1 (p) + F 2 (p), 

JO JO Jb 



hit) := 



and li 7 ^^)! < e bp 5, where \f(t)\dt < 5. Therefore, the contribution of 
the "tail" f b (t) of /, 

0, t < b, 
f(t), t>b, 

can be considered as noise if b > is large and 5 > is small. We assume 
in d2|) that / G L 2 [0,oo). One may also assume that / £ L 1 ^, 00), or that 
\f(t)\ < cie C2t , where c\,C2 are positive constants. If the last assumption 
holds, then one may define the function g(t) := /(t)e~ (c2+1)t . Then g(t) £ 
L 1 [0,oo), and its Laplace transform G(p) = F(p + C2 + 1) is known on the 
interval [02 + 1, C2 + 1 + b] of real axis if the Laplace transform F(p) of f(t) is 
known on the interval [0, b]. Therefore, our inversion methods are applicable 
to these more general classes of functions / as well. 

The operator C : Xq^ — > L 2 [0, 00) is compact. Therefore, the inversion 
of the Laplace transform ([1]) is an ill-posed problem (see [T7], (20] )■ Since the 
problem is ill-posed, a regularization method is needed to obtain a stable in- 
version of the Laplace transform. There are many methods to solve equation 
([TJ stably: variational regularization, quasisolutions, iterative regularization 
(see e.g, [13], p2], [20], [21] )• In this paper we propose an adaptive iterative 
method based on the Dynamical Systems Method (DSM) developed in [20] , 
[21j . Some methods have been developed earlier for the inversion of the 
Laplace transform (see [2], [5], [8], [E]). In many papers the data F(p) are 
assumed exact and given on the complex axis. In [16] it is shown that the 
results of the inversion of the Laplace transform from the complex axis are 
more accurate than these of the inversion of the Laplace transform from the 
real axis. The reason is the ill-posedness of the Laplace transform inversion 
from the real axis. A survey regarding the methods of the Laplace trans- 
form inversion has been given in [5]. There are several types of the Laplace 
inversion method compared in [5]. The inversion formula for the Laplace 
transform is well known: 

fit) = ^~ / F(p)J*dp, a>0, (3) 
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is used in some of these methods, and then fit) is computed by some quadra- 
ture formulas, and many of these formulas can be found in [6] and |15j . 
Moreover, the ill-posedness of the Laplace transform inversion is not dis- 
cussed in all the methods compared in [5J. The approximate f(t), obtained 
by these methods when the data are noisy, may differ significantly from f(t). 
There are some papers in which the inversion of the Laplace transform from 
the real axis was studied (see [I], [3j, [7J, [TO], [J6], [J8], [J9], [23], [21]). In 
[1] and [19] a method based on the Mellin transform is developed. In this 
method the Mellin transform of the data F(p) is calculated first and then 
inverted for f(t). In [3] a Fourier series method for the inversion of Laplace 
transform from the real axis is developed. The drawback of this method 
comes from the ill-conditioning of the discretized problem. It is shown in 
[4] that if one uses some basis functions in Xq^, the problem becomes ex- 
tremely ill-conditioned if the number m of the basis functions exceeds 20. 
In [10] a reproducing kernel method is used for the inversion of the Laplace 
transform. In the numerical experiments in [TO] the authors use double and 
multiple precision methods to obtain high accuracy inversion of the Laplace 
transform. The usage of the multiple precision increases the computation 
time significantly which is observed in [10] . so this method may be not effi- 
cient in practice. A detailed description of the multiple precision technique 
can be found in [9] and [UJ . Moreover, the Laplace transform inversion 
with perturbed data is not discussed in [10J. In [24] the authors develop 
an inversion formula, based on the eigenfunction expansion for the Laplace 
transform. The difficulties with this method are: a) the inversion formula is 
not applicable when the data are noisy, b) even for exact data the inversion 
formula is not suitable for numerical implementation. 

The Laplace transform as an operator from Cq^ into L 2 , where Cofc = 
{/(t) € C[0,+oo) | supp/ C [0,Jfe)}, k = const > 0, L 2 := L 2 [0,oo), is 
considered in [TJ. The finite difference method is used in [7J to discretize 
the problem, where the size of the linear algebraic system obtained by this 
method is fixed at each iteration, so the computation time increases if one 
uses large linear algebraic systems. The method of choosing the size of the 
linear algebraic system is not given in [7j. Moreover, the inversion of the 
Laplace transform when the data F(p) is given only on a finite interval [0, d], 
d > 0, is not discussed in [7J. 

The novel points in our paper are: 

1) the representation of the approximation solution (|73H of the function 
f(t) which depends only on the kernel of the Laplace transform, 

2) the adaptive iterative scheme (I76p and adaptive stopping rule (|87l) . 
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which generate the regularization parameter, the discrete data Fg(p) 
and the number of terms in (|73p . needed for obtaining an approxima- 
tion of the unknown function f(t). 

We study the inversion problem using the pair of spaces (Xq^, L 2 [0, d]), 
where Xq j is defined in ([2]), develop an inversion method, which can be easily 
implemented numerically, and demonstrate in the numerical experiments 
that our method yields the results comparable in accuracy with the results, 
presented in the literature, e.g., with the double precision results given in 
paper [TP] . 

The smoothness of the kernel allows one to use the compound Simp- 
son's rule in approximating the Laplace transform. Our approach yields a 
representation ()T3f) of the approximate inversion of the Laplace transform. 
The number of terms in approximation (I73D and the regularization parame- 
ter are generated automatically by the proposed adaptive iterative method. 
Our iterative method is based on the iterative method proposed in [14j . 
The adaptive stopping rule we propose here is based on the discrepancy- 
type principle, established in [22J. This stopping rule yields convergence of 
the approximation (|73p to f(t) when the noise level <5 — ► 0. 

A detailed derivation of our inversion method is given in Section 2. In 
Section 3 some results of the numerical experiments are reported. These 
results demonstrate the efficiency and stability of the proposed method. 



2 Description of the method 

Let / S Xqj,. Then equation (JT]) can be written as: 

0C/)(p) := [ b e- pt f(t)dt = F(p), < p. (4) 
J o 

Let us assume that the data F(p), the Laplace transform of /, are known 
only for < p < d < oo. Consider the mapping C m : L 2 [0, b] — ► M m+1 , where 

r b 

(An/),-:= e-Pi t f(t)dt = F(p i ), i = 0,1,2,..., m, (5) 
J o 

Pi := ih, i = 0, 1, 2, . . . , m, h := — , (6) 

m 

and m is an even number which will be chosen later. Then the unknown 
function f(t) can be obtained from a finite-dimensional operator equation 



4 



©. Let 



(u,v) W m := y]w^UjVj an d \\u\\w m '■= (u, u) 

3=0 



W< 



(7) 



be the inner product and norm in M m+1 , respectively, where Wj are the 
weights of the compound Simpson's rule (see [BJ p.58]), i.e., 



(m) 



h/3, j = 0, m; , 
Ah/3, j = 21-1, 1 = 1,2,..., m/2; h = —, (8) 
k 2h/3, j = 2l, 1 = 1,2,. ..,(m-2)/2, 



m 



where m is an even number. Then 



(£m9,v} wm =Y, w t ] / e~ Pjt 9(t)dt Vj 

3=0 J ° 

= / 9WE W ! m) e-^ t v ] dt = {g,C m v) XQb , 

J „■ n 



where 



i=o 



Vl 

\ v m J 



pm+1 



and 



and 



( Ioe-^ET=o^ m) ^ P ^dt \ 
\fSe-^ET=o^ m) ^Sdt) 



0) 



(10) 



(11) 



<9>h)xo, b : = / 9(t)h(t)dt. 

J 

It follows from ^ and CE]) that 

(C£ m5 )(t) = Y, w< P )e ~ Pjt I e- p > z 9(z)dz := (T^g)(t), (12) 



:= Q (m V (13) 
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where 

(qM) :=ti) M [ b e-te+Pi) t dt = w\ m) 1 ~ e ~ b(Pi+Pi \ i,j = 0,1,2,..., m 
' Jo J Pi+Pj 

(14) 

Lemma 2.1. Ze£ it/" 1 ^ 6e defined in ©. T/ten 



E™f ) = d > ( 15 ) 

3=0 



for any even number m. 

Proof. Prom definition ([8]) one gets 



m m/2 (m-2)/2 

E(m) (m) . (m) , \ (m) . \ -» (rr 

j=0 j=l 3=1 



2h ^4fc (m ^ /2 2/l (16) 

2/i 2/im /i(m - 2) d 

= 1 1 = rim = — m = d. 

3 3 3 m 

Lemma l2.1l is proved. □ 

Lemma 2.2. The matrix Q^ m \ defined in (I14p . is positive semidefinite and 
self-adjoint in W n+l with respect to the inner product ([7]). 

Proof. Let 

{H m )ij := / e-te+f&dt = ■ , (17) 

JO Pi + Pj 

and 

(A»)y={ ^ ' \" J; . (18) 

I U, otherwise, 

u>j m ' ) are defined in ([8j). Then (D m H m D m u,v) Rm +i = (u, D m H m D 
where 

m 

(n,w) Rm+ i := ^ujVj, u,v € M m+1 . (19) 
j=0 
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(20) 



We have 

m m 

(Q {m) u,v) wm = Y,w { 3 m) (Q^u) 3 v 3 = Y,(DmH m D mU ) ]Vj 

3=0 j=0 

= {D m H m D m u, f)jj m +i = (u, D m H m D m v) jgrn+i 

m m 

= ^ Uj (D m H m D m v) j = y~] Uj wf 1 ' (H m D m v) j 

j=0 j=0 

= {u,Q^v) Wm . 

Thus, is self-adjoint with respect to inner product ([7]). We have 

(H m )ij= I* e-te+P&dt = t e- pit e- p ^dt 

K >3 Jo Jo (21) 

where {'r)x 0b ^ s defined in (jlip . This shows that H m is a Gram matrix. 
Therefore, 

(H m u,u) Rm+ i > 0, \/u G R m+1 . (22) 

This implies 

(Q( m >u,u) wm = (Q (m) n,D m n) Rm+ i = (H m D m u, D m u) Rm+1 > 0. (23) 

Thus, Q( m ) is a positive semidefinite and self-adjoint matrix with respect to 
the inner product ([7]). □ 

Lemma 2.3. Let be defined in (|12p . Then T^ m ' is self-adjoint and 

positive semidefinite operator in Xq^ with respect to inner product (llip . 

Proof. From definition (|12p and inner product (lllh we get 

.fern „{, 
(TWg,h) Xob = Y, w P e ~ Pjt e~K*g{z)dzh(t)dt 

J j_Q 

= f 9(*)Y,Wj n) e~ Pi ' I e~^ l h{t)dtdz (24) 



3=0 

(g,T^h) 



Xn 



Thus, r( m ) is a self-adjoint operator with respect to inner product (jlip . Let 
us prove that T( m ) is positive semidefinite. Using jT2]), ©, © and (dU, 



one gets 



(T^g,g) 



Xn 



/ y"^ m) e" Pj * / e- p > z g(z)dzg(t)dt 
Jo j=0 Jo 

m b b 

Viflf 1 / e- p i z g{z)dz / e- p ^g(t)dt 
~Z Jo Jo 



(25) 



3=0 
in 



j=0 KJ0 



e~^ z g{z)dz) > 0. 



□ 



Lemma 12.31 is proved. 

From (fTUj) we get Range[£^] = spanlvjj 71 ^ k(pj, •, 0)}jL , where 

k(p,t,z) := e - p{t+z \ 
Let us approximate the unknown f{t) as follows: 

/(«) * £cf = T~iC m F {m) := / m (t), 

3=0 

where pj are defined in ([6]), T aim is defined in (f34"|) . and cj- are constants 
obtained by solving the linear algebraic system: 



(26) 



(27) 



(a/ + Q (m) )c (m) =F ( - m \ 
where Q( m ) is defined in (|13p . 



(28) 



,(m) ._ 



and F< m > := 



/ F(po) \ 
F{ P i) 

\ F( Pm ) J 



(29) 



To prove the convergence of the approximate solution f(t), we use the fol- 
lowing estimates, which are proved in |21| . so their proofs are omitted. 

Lemma 2.4. Let and be defined in (|12|) and (|13|) . respectively. 

Then, for a > 0, £/ie following estimates hold: 



\\Qa,in^m\\ < 



2^0"' 



(30) 
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4Qain\\ < h (31) 

\\T~L\\ < -, (32) 
a 

ll-^ajm^mll — (^) 

where 

Q a , m ■= Q {m) + al T a , m :=T^+aI, (34) 

/ is the identity operator and a = const > 0. 

Estimates (f30|) and (f3Tj) are used in proving inequality (f92|) . while esti- 
mates ([32]) and (j33|) are used in the proof of lemmas 2.9 and 2.10, respec- 
tively. 

Let us formulate an iterative method for obtaining the approximation 
solution of f(t) with the exact data F(p). Consider the following iterative 
scheme 

u n {t) = gu n _i(t) + (1 - q)T££*F, u (t) = 0, (35) 
where C* is the adjoint of the operator C, i.e., 

(C*g)(t)= [ d e-P t g(p)dp, (36) 



(37) 



(Tf)(t) := (£*£/)(*) = / [ k(p,t,z)dpf(z)dz 



o jo 



/(*) 



t + Z 

k(p,t,z) is defined in ([26]) . 

T a :=al + T, a > 0, (38) 

an := gan-i, a > 0, g G (0, 1). (39) 

Lemma 2.5. Xei T a be defined in (j38|) . £/ = F, and / _L N{C), where 
M(C) is the null space of C Then 

a\\T- l f\\ as a 0. (40) 

Proof. Since / _L J\f(£), it follows from the spectral theorem that 

f°° a 2 

hma 2 ||T-7|| 2 = lim / -——d(E s f,f) = \\P mc) f\\ 2 = 0, 

where E s is the resolution of the identity corresponding to £*£, and P is 
the orthogonal projector onto M{C). 

Lemma 12.51 is proved. □ 
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Theorem 2.6. Let Cf = F, and u n be defined in (|35p Then 

lim ||/-un|| =0. (41) 

n—>oo 

Proof. By induction we get 

n-l 

u n = J2^ n % 1 +1 ^F, (42) 

3=0 

where T a is defined in (|55|) . and 

(n) . = g n-j-l _ q n-j ^ 



3 

Using the identities 



Cf = F, (44) 
T- X CC = T- 1 {T + aI-aI) = I- aT' 1 (45) 



and 



we get 



n-l 

,(») 



5>} n) = l-? n , (46) 

3=0 



Therefore, 



n—l n— 1 

3=0 i=0 

n-l 

3=0 



71-1 

■ (»), 



(47) 



11/ - u n \\ < q n \\f\\ + £>f Vlll^/ll- (48) 

3=0 

To prove relation (|4ip the following lemma is needed: 

Lemma 2.7. Let g(x) be a continuous function on (0, oo), c > and q G 
(0, 1) 5e constants. If 

\im g(x) = g(0) := g , (49) 

n-l 

lim V (q^- 1 - q n -i) g{cq^ +1 ) = g . (50) 

3=0 
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Proof. Let 



1-1 



(51) 



where wj n ' ) are defined in (|43p . Then 



|F n +i(n)-5o| < l^/WI + 



s;«i B W +i )-fld 



Take e > arbitrarily small. For sufficiently large fixed /(e) one can choose 
n(e) > /(e), such that 

|*i( e )(n)|<|, Vn>n(e), 

because lim^^oo q n = 0. Fix I = 1(e) such that \g(cq J ) — go I < § for j > /(e). 
This is possible because of dM|). One has 



\F l[e) (n)\ <-,n> n(e) > /(e) 



and 



£ ^fW +1 )-90 



< £ u,fW +1 )-<?o| + | £ ^-%o| 



(n) 



J=*(e) 



i=*(e) 



4 

< I + M<Z n ~' (e) < e, 
if n(e) is sufficiently large. Here we have used the relation 



l>f = i-<r-*. 

Since e > is arbitrarily small, relation (I50p follows. 
Lemma 12.71 is proved. 

Lemma 12.51 together with Lemma 12.71 with g(a) = aUT" 1 /!! yield 



n-l 



feE w i B Wii^/ii = o. 



□ 



(52) 



3=0 
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This together with estimate ([4*8]) and condition q £ (0, 1) yield relation (|4ip. 
Theorem 12.61 is proved. □ 

Lemma 2.8. LetT andT^ be defined in (I37p and fll2|> . respectively. Then 

{2bdf 



\\rp _ ji(m)|| < 

~ 540\/T0m 4 " 
Proof. From definitions (|37p and (|12p we get 



(53) 



|(T-TM)/(i)|< / / fcfot^dp-V^H^M) 



\f{z)\dz 



< 



180m 4 P e[o,d] 
b d 5 



max (i + z) e 



4 -p(t+z) 



o 180m 4 



(t + zf\f(z)\dz< 



d 5 



180m 4 



(t + 6) 9 - t 9 
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1/2 



180m 4 



A,, 



1/2 



(54) 



where the following upper bound for the error of the compound Simpson's 
rule was used (see [6J p. 58]): for / £ C^[xq,X2z], a?o < x 2 i, 



f X21 h 
/ f(x)dx - - 



I l-l 
/o+4^/ 20 _ 1) + 2^/ 2i + / : 
i=i i=i 



J-21 



< Rh (55) 



where 



fj '■= f( x j)> x j = xo+jh, j = 0,1,2,. . . ,2/, /i 



X2l ~ x 
21 



and 



Ri 



(X 2 l - Xpf 

180(2/) 4 



l/ (4) (£)l, x <i<x 2l . 



(56) 
(57) 



This implies 



\\{T-T^)f\\ Xofi < 



d r > 



540m 4 



(2b) 10 - 2fc 10 
10 



1/2 



An h < 



(2bd) 1 



540v / 10m 4 



so estimate (j53j) is obtained. 
Lemma 12.81 is proved. 



Ao,b j 
(58) 

□ 
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Lemma 2.9. Let < a < qq, 



m = k 



K>0. 



Then 



irp _ rp{rn)\\ < 



(2bd) 1 



(59) 
(60) 



540VT0a K 4 

where T and are defined in (13TI) and (j!2j) . respectively. 
Proof. Inequality (|60p follows from estimate (153p and formula (I59h . □ 
Lemma 12.91 leads to an adaptive iterative scheme: 

U n ,m n (t) = qUn-^m,,-! + (1 ~ ^^m™^^"^ ' «0,mo(*) = °> ( 61 ) 



where g € (0,1), a n are defined in (j39j) . T 0jm is defined in (f34l) . vl m £ is 
defined in (|5|), and 



( F(po) \ 
F(pi) 

\ F(Pm) J 



om+1 



(62) 



Pj are defined in (|6|). In the iterative scheme (|6ip we have used the finite- 
dimensional operator approximating the operator T. Convergence of 
the iterative scheme ([6Tj) to the solution / of the equation Cf = F is estab- 
lished in the following lemma: 

Lemma 2.10. Let Cf = F and u n , mn be defined in (foTj) . If m n are chosen 
by the rule 



1/4 



a n = qctn-l, 9 6(0,1), K,a >0, (63) 



where \[x]~\ is the smallest even number not less than x, then 

lim \\f-u n , mn \\ =0. (64) 

n^oo 

Proof. Consider the estimate 



11/ u n,m n || — 



Un\\ + «n - U 



n u, n,m n 



:= h(n) + I 2 (n), (65) 



13 



where I\{n) := \\f — u n \\ and h^n) := \\u n — ii n ,mj|- By Theorem 12.61 
we get I\{n) — > as n — > oo. Let us prove that lim n _>oo = 0. Let 
u n iin,m n - Then, from definitions f|35 j) and (|6ip . we get 

U n = qU n . 1 + (l-q)(T- n 1 C*F-T-] mn C* mn F( m ^, U = 0. (66) 

By induction we obtain 

n-l 

I7 n = J>f {T'l^F - T-l umj+i (C mj+1 TF^) , (67) 

3=0 

where Uj are defined in ([4*3j) . Using the identities £/ = F, £ m / = F^ m \ 

T- X T = T-\T + a I-aI)=I- aT' 1 , (68) 
T-i TM = T-^(T( m ) + a/ - a/) = / - or-i , (69) 

t-1 - r- 1 = r-^(r - tHjt; 1 , (70) 



one gets 



n-l 
n-l 

= £a,f a i+1 r-| limj . +l ( T - T(mj+1) ) T 4,/- 

J'=0 



(71) 



This together with the rule (|63p . estimate (|32p and Lemma 12.81 yield 



n-l 

,(n) 



lic/nll < ^II^^Jir-T^+^iliir-^/ii 



(72) 



(2bdf ^ {n) 



< 



J2u^a j+1 \\T-lJ\\. 



540A/l0a K 4 j =Q 

Applying Lemma 12.51 and Lemma 12.71 with g(a) = a\\T~ f\\, we obtain 
lim^oo \\U n \ \ = 0. 

Lemma 12.101 is proved. □ 
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2.1 Noisy data 

When the data F(p) are noisy, the approximate solution ([77]) is written as 



3=0 
Jm,S) 



rp-l f* pW 

a,m J ~'m 5 ' 



(73) 



where the coefficients cj- ' are obtained by solving the following linear 
algebraic system: 

Q a m is defined in 



(m,S) 



7 (m) 
6 ' 



(74) 



i c \ 

(m,<5) 



F 



(m) ._ 



/ *«(po) \ 
^(pi) 



(75) 



w 



(m) 



are defined in ([8]) , and pj are defined in 



To get the approximation solution of the function f(t) with the noisy 
data F$(p), we consider the following iterative scheme: 



? U tl,ra„_i + (1 - (l) T a^m n £m n F c 



0, 



(76) 



S ' u 0,m 

where T am is defined in (|34"|) . a n are defined in ([39]) . q 6 (0,1), -^i is 
defined in ([75]) , and m n are chosen by the rule ([63]) . Let us assume that 

Fs( Pj )=F(p J ) + 6 j , 0<|<5 i |<<y, j = 0,l,2,...,m, (77) 

where Sj are random quantities generated from some statistical distributions, 
e.g., the uniform distribution on the interval [—5,5], and 5 is the noise level 
of the data F(p). It follows from assumption (|77p . definition ([8]), Lemma [2.1l 
and the inner product ([7]) that 

m m 
3=0 3=0 

Lemma 2.11. Letu n ^ mn andu 5 nmn be defined in ([6T]) and (f76|) . respectively. 
Then _ 

X <l6 -.(l-q n ), 9 €(0,1), (79) 



*n,m„ u n,m n II — 2 



where a n are defined in 
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Proof. Let U* := u n,m n ~ u nm n - Then, from definitions (|6ip and (|76p . 

Ui = qUi_ l + {l- q )T^ mn L* m SF^-F^ n \ U* = 0. (80) 
By induction we obtain 



n-l 

U 5 n =J2^ T ati,rn j+l (£rn j+1 T(F^ - i^ +l) ), (81) 
j=0 



where tu^ are defined in fi3]) . Using estimates ([78]) and inequality ([33]) . 
one gets 

where are defined in ([33"]) . 

Lemma 12.111 is proved. □ 

Theorem 2.12. Suppose that conditions of Lemma \2.10\ hold, and ns sat- 
isfies the following conditions: 

lim nx = oo, lim — ■ = 0. (83) 

Then 

lim||/-< )mn J| = 0. (84) 
Proof. Consider the estimate: 

11/ - u i 6 , mns II < 11/ - u n s ,m ns II + \\Un s ,mn s ~ ^n,,m„, II- ( 85 ) 
This together with Lemma 12.111 yield 

11/ " < || < ||/ - u ne>mn6 1| + -^(1 - q n ). (86) 

Applying relations ([83]) in estimate ([86]) , one gets relation ([51]) . 

Theorem 12.121 is proved. □ 

In the following subsection we propose a stopping rule which implies 
relations ([83]) . 
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2.2 Stopping rule 

In this subsection a stopping rule which yields relations (|83|) in Theorem l2.12l 
is given. We propose the stopping rule 

G ns , mns < C5 £ < G n , mn , 1 < n < n s , C > Vd, e £ (0, 1), (87) 

where 

G n , mn = qG n ^ mn _ x + (1 - q)\\C m J m ^ - Fl mn) \\ wmn , G 0imo = 0, (88) 
|| • \\w m is defined in ©, 

m 

.M^^^^Jg-Pit (89) 

3=0 

Wj and pj are defined in ([8]) and (J6]) , respectively, and c^ 1 ^ are obtained 
by solving linear algebraic system (|74p . 
We observe that 

/- ~(mn,5) _ p( m ") _ n( m ")^( m «- 5 ) _ p( m ") 
t -ra n z r g — V c ^5 

= Q^(a n I + Ql" 1 "))- 1 ^™"' - 

= (Q {mn) + a n I - a n I){a n I + Q(™n))-i F( (™») - 

= -a n (a n I + Q^) r i F {m n ) = _ anC (m n ,5)^ 

(90) 

Thus, the sequence (|88p can be written in the following form 

Gn,m n = qGn-l,^.! + (1 ~ ?)o„ || C (m " ' 5) || f™. , G ,m = 0, (91) 

where || • \\w m is defined in ([7]), and c^ m ' 5 ^ solves the linear algebraic systems 
(EH). 

It follows from estimates (1751). (1301) and (I3T1) that 



fln||c (m "' 5) |k™n =a n ||(a n / + Q( m "))- 1 F ( { m " ) || wm n 

< anlK^i + Q^))- 1 ^^-^)!!^ 

+ a n ||(a n / + Q( m "))- 1 F( m ")|| H , m „ 

< [| F j m «)_ F (™«)|| wmn 

+ a n ||(a n / + Q( m «))- 1 £ mn /||^ m „ 
<5>/d + v^[|/l|jfo6- 



(92) 
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This together with (|9T|) yield 

+ 0—q){sVd + ^:\\f\\ XOib ), (93) 

or 

G n , mn -5Vd< q(G n -i :rrin _ 1 - 5Vd) + (1 - q)V^\\f\\ Xo>b . (94) 
Lemma 2.13. The sequence (|9ip satisfies the following estimate: 



G„,„,,-^< (1 -^" /lk " , (95) 



where a n are defined in (|39j) . 
Proof. Define 

:= G n , mn - <fVd (96) 

and 

^n:=(l-?)v^ll/lk,»- (97) 
Then estimate (|94p can be rewritten as 

*n<9*n-l + v^n-l, (98) 

where the relation a n = qa n -\ was used. Let us prove estimate ([95]) by- 
induction. For n = we get 

* Q = -5Vd< {1 - 9) ^ mx <>\ (99) 
1 - v 7 ^ 

Suppose estimate (i95j) is true for < n < k. Then 

q 

^k+i < q^k + Vqipk < -, -pipk + Vqipk 

i - V? 

^ = i-.^„ ^i (ioo) 



i - y/q i - \fq ^k+i 

Vq yfOk , i , 

-_ipk+i = - 



1 - y/q V«fc+i 1 - V? 

where the relation a^+i = qa^ was used. 

Lemma 12.131 is proved. □ 
Lemma 2.14. Suppose 

Gi, mi > dVd, (101) 



where G njfnn are defined in (|9T|) . T/ien i/iere exis£ a unique integer n$, 
satisfying the stopping rule (|87|) wii/i C > \/cZ. 
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Proof. Prom Lemma 12.131 we get the estimate 



G n ,m n < 6Vd-\ -, (102) 

1 - y/q 



where a n are defined in (|39|) . Therefore, 

limsup G ntmn < 5 Vd, (103) 

n^oo 

where the relation lim^^co a n = was used. This together with condition 
(jlOip yield the existence of the integer n$ . The uniqueness of the integer ns 
follows from its definition. 

Lemma 12.141 is proved. □ 

Lemma 2.15. Suppose conditions of Lemma \2. 1J\ hold and ns is chosen by 
the rule (I87p . Then 

lim = 0. (104) 

Proof. From the stopping rule (|87p and estimate f)102|) we get 

C5 <G ns - 1 . mn/> _ 1 <5Vd+ ' ' > 

where C > yd, e G (0, 1). This implies 

SjCS*- 1 - y/d) < (l-q)\\f\\x , b 

so, for e G (0, 1), and a ns = qa ns -i, one gets 

lim — ■ = lim ^ < lim — — jj^jj^t _ q (107) 

S^o ja^ S^o _ q )(C - 5 l - £ \/d) 

Lemma l2.15l is proved. □ 

Lemma 2.16. Consider the stopping rule (|87p . where the parameters m n 
are chosen by rule (f63j) . If ns is chosen by the rule (|87l) then 

lim rii = oo. (108) 

<5^o 
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Proof. From the stopping rule flSTJ) with the sequence G n defined in (f9T|) 
one gets 

qC 5 e + (1 - 9)0^11^^11^ < 9Gn J -l,m n5 _ 1 
+ (1 - 9KJc( m ^|| w -n a = G n , < C8 e , 



where c^ m '^ is obtained by solving linear algebraic system (|74p . This implies 

0<a n4 ||c( m »*»^|| wm n S <C<f. (110) 

Thus, 

lima n4 || C ( m ^^||^ =0. (Ill) 
0— »o 

If / 0, then there exists a > such that 

^0 

where £^ m ' 1 is the resolution of the identity corresponding to the operator 
Q (m) := Let 

MM := a 2 ||g-^i m) || 2 vi, m , Q m ,a := al + Q™. 
For a fixed number a > we obtain 



h m (5,a)=a 2 \\Q-] a Ft ) » 2 



00 a 2 



> 



(a + s) 2 
A <> a 2 



a 

> 



(« + *) 2 ^ " ' ° ' (113) 

(m) 



2 „\( m ) 



(a + A ) 2 7o 

-2||p( m ) 77( m )||2 

(« + Ai m) ) 2 

Since E^' is a continuous operator, and ||i^( m ) — F^ \\w m < VdS, it follows 
from (fTT2j) that 

Em(£j)^),fjV = (^F^^M)^ > 0. (114) 
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Therefore, for the fixed number a > we get 

h m (5,a) > c 2 > 



(115) 



for all sufficiently small 5 > 0, where c 2 is a constant which does not depend 
on 5. Suppose lim^o a n& ^ 0. Then there exists a subsequence Sj — » as 
j — > oo, such that 



and 
< m 



a ns . > ci > 0, 



^(aoK,,) 1/4 ]l < [Kao/ci) 1/4 ] 



(116) 



:= C3 < oo, k, ao > 0, 
(117) 

where the rule f)63|) was used to obtain the parameters m ns . This together 
with (HHP and (fTToll yield 



(m„ 5 .) K { ) 



lim /i m (<$,-, a„ 5 .) > lim 



"na . ii"" (win. . ) 8j 

An J 



Kj ) 

Or- + A J ) 2 



J||F (W V p( m ^)||2 

> liminf °- ; ; > 0. 



(118) 



K. ) 



This contradicts relation (jllip . Thus, lim^o a n a = lim^o a-oQ™ 6 = 0, i.e., 
lim^o ns = oo. 

Lemma l2.16l is proved. □ 



It follows from Lemma 12.151 and Lemma 12.161 that the stopping rule (|87h 
yields the relations (|83p . We have proved the following theorem: 



Theorem 2.17. Suppose all the assumptions of Theorem \2.12\ hold. m n are 
chosen by the rule (|63|) . ng is chosen by the rule (|87l) and Gi, mi > C5, where 
G n ,m n are defined in (f9Tj) . then 



limJI/-< 



5^0 



6,™,ng I 



0. 



(119) 



2.3 The algorithm 

Let us formulate the algorithm for obtaining the approximate solution f^: 
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(1) The data Fg(p) on the interval [0, d], d > 0, the support of the function 
f(t), and the noise level 5; 

(2) initialization : choose the parameters k > 0, ao > 0, q G (0,1), e 6 
(0, 1), C > and set u 5 Qmo = 0, G Q = 0, n = 1; 

(3) iterate, starting with n = 1, and stop when condition (|126h ( see below) 
holds, 

(a) a n = a </\ 

(b) choose m n by the rule (|63l) , 

(c) construct the vector F^ 1 "" 1 : 

(pjmn)^ = Fs{pi) ^ ^ = lK h = d/m ^ / = o, 1, . . . ,m, (120) 

(d) construct the matrices H mn and D mn : 

{H mn ) ir .= e -^ +p ^dt = ■ , i,j = 1,2,3,... ,m n 

Jo Pi+ Pj 

(121) 

( (m n ) ■ _ . 

m?l lJ I 0, otherwise, 

where u?j m ^ are defined in ©, 

(e) solve the following linear algebraic systems: 

(a n J + H mn D mn )c^ m ^ = (123) 

where = Cj (mn ' 5) , 

(f) update the coefficient c j m "' <5 ' 1 of the approximate solution u\ i mn (t) 
defined in (|73p by the iterative formula: 



< mn (t) = ^_ 1 , mn _ 1 (i) + (1 - ff^c^^e-*', (124) 



J=l 

where 

4, mo (t) = 0- (125) 
Stop when for the first time the inequality 

G n , mn = qG n - lt m n ^+a n \\c^ s \\ Wmn < C5 £ (126) 

holds, and get the approximation f s (t) = u ^ (t) of the function 
/(*) by formula ([124]l . 
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3 Numerical experiments 



3.1 The parameters k, ao, d 

From definition ([39]) and the rule (f63j) we conclude that m n — > oo as a n — > 0. 
Therefore, one needs to control the value of the parameter m n so that it 
will not grow too fast as a n decreases. The role of the parameter k in 
(|63p is to control the value of the parameter m n so that the value of the 
parameter m n will not be too large. Since for sufficiently small noise level 5, 
namely 5 £ (10 -16 , 10~ 6 ], the regularization parameter a ng , obtained by the 
stopping rule (|87|) . is at most O(10~ 9 ), we suggest to choose k in the interval 
(0,1]. For the noise level 5 E (10~ 6 ,10~ 2 ] one can choose k E (1,3]. To 
reduce the number of iterations we suggest to choose the geometric sequence 
a n = ao5 an , where ao £ [0.1, 0.2] and a £ [0.5, 0.9]. One may assume without 
loss of generality that 6=1, because a scaling transformation reduces the 
integral over (0, 6) to the integral over (0, 1). We have assumed that the data 
F{p) are defined on the interval J := [0, d]. In the case the interval J = 
[di,d], < d\ < d, the constant d in estimates JBDJ), (JTSJ) , ((79j) , d82j), AMI) . 
([95]) . and (|102p are replaced with the constant d — d\. If 6 = 1, i.e., f(t) = 
for i > 1, then one has to take d not too large. Indeed, if /(£) = for £ > 1, 
then an integration by parts yields: F(p) = [/(0) — e~ p f(l)]/p + 0(l/p 2 ), 
p — > oo. If the data are noisy, and the noise level is 5, then the data becomes 
indistinguishable from noise for p = 0(1/5). Therefore it is useless to keep 
the data F$(p) for d > 0(1/ 5). In practice one may get a satisfactory 
accuracy of inversion by the method, proposed in this paper, when one uses 
the data with d £ [1, 20] when 5 < 10~ 2 . In all the numerical examples we 
have used d = 5. Given the interval [0,d], the proposed method generates 
automatically the discrete data F$(pj), j = 0, 1, 2, . . . , m, over the interval 
[0, d] which are needed to get the approximation of the function f(t). 

3.2 Experiments 

To test the proposed method we consider some examples proposed in [I], 
0, [2], [4], [5], [8], [ID], pi], [18] and [24] ■ To illustrate the numerical 
stability of the proposed method with respect to the noise, we use the noisy 
data F$(p) with various noise levels 5 = 10~ 2 , 5 = 10~ 4 and 5 = 10 -6 . 
The random quantities 6j in (|77p are obtained from the uniform probability 
density function over the interval [— S, 5]. In examples 1-12 we choose the 
value of the parameters as follows: a n = 0.1q n , q = and d = 5. The 
parameter k = 1 is used for the noise levels 5 = 10~ 2 and 5 = 10~ 4 . When 
5 = 10 -6 we choose k = 0.3 so that the value of the parameters m n are not 
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very large, namely m n < 300. Therefore, the computation time for solving 
linear algebraic system (|123p can be reduced significantly. We assume that 
the support of the function f(t) is in the interval [0, b] with b = 10. In the 
stopping rule ([57)1 the following parameters are used: C = \fd + 0.01, e = 
0.99. In example 13 the function f(t) = e~ t is used to test the applicability 
of the proposed method to functions without compact support. The results 
are given in Table 13 and Figure 13. 

For a comparison with the exact solutions we use the mean absolute 
error: 



MAE :-- 



E?=i(f(ti)-fL 5 (ti)f 
100 



1/2 

, ^ = 0.01+0.1(3-1), j = l,...,100, 



(127) 

where f(t) is the exact solution and (t) is the approximate solution. 
The computation time (CPU time) for obtaining the approximation of fit), 
the number of iterations (Iter.), and the parameters m ns and a ns generated 
by the proposed method are given in each experiment (see Tables 1-12). All 
the calculations are done in double precision generated by MATLAB. 



Example 1. (see [10]) 

f 1, 1/2 < * < 3/2, „ , , f 1, P = 0, 

/l(t) = \0, otherwise, Fl{p) = p>0 . 




Figure 1: Example 1: the stability of the approximate solution 



The reconstruction of the exact solution for different values of the 
noise level 5 is shown in Figure 1. When the noise level 5 = 10~ 6 , our 
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Table 1: Example 1. 



5 




MAE 




m ns 


Iter. 


CPU time (second) 






1.00 x 10" 


-2 


9.62 x 10" 


2 


30 


3 


3.13 x 10~ 2 


2.00 x 10" 


-3 


1.00 x 10" 


-4 


5.99 x 10" 


2 


32 


4 


6.25 x 10~ 2 


2.00 x 10" 


-7 


1.00 x 10" 


-6 


4.74 x 10- 


2 


54 


5 


3.28 x 1Q- 1 


2.00 x 10" 


-10 



result is comparable with the double precision results shown in [10J. 
The proposed method is stable with respect to the noise 5 as shown 
in Table 1. 

• Example 2. (see J777|/ ) 




p = 0, 
p > 0. 



2 4 6 8 10 



Figure 2: Example 2: the stability of the approximate solution 



Table 2: Example 2. 



5 


MAE 




Iter. 


CPU time (seconds) 


Q"ng 




1.00 x lO" 2 


1.09 x lO" 1 


30 


2 


3.13 x 10~ 2 


2.00 x 10" 


-3 


1.00 x 10~ 4 


8.47 x 10~ 2 


32 


3 


6.25 x 10~ 2 


2.00 x 10" 


-6 


1.00 x 10" 6 


7.41 x 10~ 2 


54 


5 


4.38 x lO" 1 


2.00 x 10" 


12 



The reconstruction of the function /2(f) is plotted in Figure 2. In [10] 
a high accuracy result is given by means of the multiple precision. But, 
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as reported in |10j . to get such high accuracy results, it takes 7 hours. 
From Table 2 and Figure 2 we can see that the proposed method yields 
stable solution with respect to the noise level 5. The reconstruction of 
the exact solution obtained by the proposed method is better than the 
reconstruction shown in [4J. The result is comparable with the double 
precision results given in [10] . For 5 = 10~ 6 and k = 0.3 the value of 
the parameter m ns is bounded by the constant 54. 

• Example 3. (see JT/, gy, J5U, /7gj/ . J&fl /j 




Figure 3: Example 3: the stability of the approximate solution 



Table 3: Example 3. 



5 


MAE 




Iter. 


CPU time (seconds) 


a ns 


1.00 x 10" 2 


2.42 x 10" 2 


30 


2 


3.13 x 10" 2 


2.00 x 10" 3 


1.00 x 10~ 4 


1.08 x 10" 3 


30 


3 


3.13 x 10" 2 


2.00 x 10~ 6 


1.00 x 10" 6 


4.02 x 10" 4 


30 


4 


4.69 x 10" 2 


2.00 x 10" 9 



We get an excellent agreement between the approximate solution and 
the exact solution when the noise level S = 10~ 4 and 10 -6 as shown in 
Figure 3. The results obtained by the proposed method are better than 
the results given in [3] . The mean absolute error MAE decreases as the 
noise level decreases which shows the stability of the proposed method. 
Our results are more stable with respect to the noise 5 than the results 
presented in [24], The value of the parameter m ng is bounded by the 
constant 30 when the noise level 5 = 10 -6 and k = 0.3. 
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• Example 4- (see 



f 1 _ o < t < 10, 

| 0, elsewhere. 




Figure 4: Example 4: the stability of the approximate solution 

As in our example 3 when the noise 5 = 10" 4 and 10~ 6 are used, we 
get a satisfactory agreement between the approximate solution and the 
exact solution. Table 4 gives the results of the stability of the proposed 
method with respect to the noise level 5. Moreover, the reconstruction 
of the function fn{t) obtained by the proposed method is better than 
the reconstruction of fA{t) shown in [5], and is comparable with the 
double precision reconstruction obtained in |10j . 



Table 4: Example 4. 



s 


MAE 




Iter. 


CPU time (seconds) 






1.00 x 1Q- 2 


1.59 x 10~ 2 


30 


2 


3.13 x 10~ 2 


2.00 x 10" 


-a 


1.00 x 10~ 4 


8.26 x 10~ 4 


30 


3 


9.400 x 10~ 2 


2.00 x 10" 


-6 


1.00 x 10" 6 


1.24 x 10" 4 


30 


4 


1.250 x 10" 1 


2.00 x 10" 


-9 



In this example when 5 = 10 6 and k = 0.3 the value of the parameter 
m ng is bounded by the constant 109 as shown in Table 4. 

• Example 5. (see g]/, \Sj) 
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f 5 (t) = 2/V3e- t/2 sm(tV3/2) 

1 - cos(10^/2)e- 10 ( p+0 - 5 > _ 2{p + 0.5)e- 1Q ( p + a5 ) sin(10vg/2) 
5(P) ~ [(p + 0.5) 2 + 3/4] " V3[(p + 0.5) 2 + 3/4] 




Figure 5: Example 5: the stability of the approximate solution 



Table 5: Example 5. 



5 


MAE 




Iter. 


CPU time (seconds) 






1.00 x 10" 2 


4.26 x 10" 2 


30 


3 


6.300 x 10~ 2 


2.00 x ltr 


-3 


1.00 x 10~ 4 


1.25 x 10~ 2 


30 


3 


9.38 x 10~ 2 


2.00 x 10" 


-6 


1.00 x 10~ 6 


1.86 x 10~ 3 


54 


4 


3.13 x 10~ 2 


2.00 x 10" 


-9 



This is an example of the damped sine function. In [2] and [8] the 
knowledge of the exact data F(p) in the complex plane is required 
to get the approximate solution. Here we only use the knowledge 
of the discrete perturbed data Fg(j)j), j = 0,1,2, ... ,m, and get a 
satisfactory result which is comparable with the results given in [2] 
and [8] when the level noise 5 = 10 -6 . The reconstruction of the 
exact solution fs(t) obtained by our method is better than this of the 
method given in [4j. Moreover, our method yields stable solution with 
respect to the noise level 5 as shown in Figure 5 and Table 5 show. In 
this example when k = 0.3 the value of the parameter m ng is bounded 
by 54 for the noise level 5 = 10" 6 (see Table 5). 

• Example 6. (see [10]) 
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' t, < t < 1, 

< 3/2 -t/2, 1 < t < 3, 
0, elsewhere. 

/ 3/2, p = 0, 

S l-e-P(l+p) | e-^+e 2 n2p-l) . n 

[ — 1 w ' p > u ' 



0.8 

0.4 
0.2 

0^ 



Figure 6: Example 6: the stability of the approximate solution 



Table 6: Example 6. 



<5 


MAE 


m ns 


Iter. 


CPU time (seconds) 




1.00 x 10" 2 


4.19 x 10" 2 


30 


2 


4.700 x 10" 2 


2.00 x 10" 3 


1.00 x 10" 4 


1.64 x 10" 2 


32 


3 


9.38 x 10~ 2 


2.00 x 10~ 6 


1.00 x 10" 6 


1.22 x 10~ 2 


54 


4 


3.13 x 10" 2 


2.00 x 10" 9 



h(t) = 

F 6 (p) = 




Example 6 represents a class of piecewise continuous functions. ^From 
Figure 6 the value of the exact solution at the points where the function 
is not differentiable can not be well approximated for the given levels 
of noise by the proposed method. When the noise level 5 = 10 -6 , our 
result is comparable with the results given in [10]. Table 6 reports 
the stability of the proposed method with respect to the noise 5. It 
is shown in Table 6 that the value of the parameter m generated by 
the proposed adaptive stopping rule is bounded by the constant 54 for 
the noise level 5 = 10 -6 and k = 0.3 which gives a relatively small 
computation time. 

• Example 7. (see [10]) 
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fr(t) = I 



Fi{p) 



-teT % - e - * + 1, < t < 1, 
l-2e -1 , 1 < t < 10, 

0, elsewhere, 



3/e- 1 + 9(1 -2/e), 

r -l-p e 1+p -e(l+P) 2 +P(3+2p) 
p(P+T? 



p = 0, 




Figure 7: Example 7: the stability of the approximate solution 



Table 7: Example 7. 



5 


MAE 




Iter. 


CPU time (seconds) 


Q"ng 


1.00 x 10" 2 


1.52 x 10" 2 


30 


2 


4.600 x 10" 2 


2.00 x 10" 3 


1.00 x 10~ 4 


2.60 x 10~ 3 


30 


3 


9.38 x 10~ 2 


2.00 x 10~ 6 


1.00 x 10" 6 


2.02 x 10~ 3 


30 


4 


3.13 x 10~ 2 


2.00 x 10~ 9 



When the noise level 5 = 10~ 4 and 5 = 10 -6 , we get numerical results 
which are comparable with the double precision results given in [10] . 
Figure 7 and Table 7 show the stability of the proposed method for 
decreasing 5. 



Example 8. (see ffl) 

fs(t) ~- 



At 2 e~ 2t , 
0, 



< t < 10, 
elsewhere. 



F S (p) 



8 + 4e~ 1 Q( 2 +p) [-2 - 20(2 + p) - 100(2 - p) 2 } 
(2+p) 3 ' 



30 




Figure 8: Example 8: the stability of the approximate solution 

The results of this example are similar to the results of Example 3. The 
exact solution can be well reconstructed by the approximate solution 
obtained by our method at the levels noise 5 = 10~ 4 and 5 = 10~ 6 
(see Figure 8). Table 8 shows that the MAE decreases as the noise 
level decreases which shows the stability of the proposed method with 
respect to the noise. In all the levels of noise 5 the computation time 
of the proposed method in obtaining the approximate solution are 
relatively small. We get better reconstruction results than the results 
shown in [JJ. Our results are comparable with the results given in [3]. 



Table 8: Example 8. 



5 MAE m ns 


Iter. 


CPU time (seconds) 




1.00 x 10" 2 2.74 x 10~ 2 30 


2 


1.100 x 10~ 2 


2.00 x 10" 3 


1.00 x 10~ 4 3.58 x 10~ 3 30 


3 


3.13 x 10~ 2 


2.00 x 10" 6 


1.00 x 10~ 6 5.04 x 10~ 4 30 


4 


4.69 x 10~ 2 


2.00 x 10~ 9 


• Example 9. (see [18]) 








h{t) = < 


' 5-t, 


< t < 5, 
elsewhere, 




*b(p) = < 


f 25/2, p = 0, 
{ s=^tl, P >0. 





As in Example 6 the error of the approximate solution at the point 
where the function is not differentiable dominates the error of the 
approximation. The reconstruction of the exact solution can be seen 
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Figure 9: Example 9: the stability of the approximate solution 

in Figure 9. The detailed results are presented in Table 9. When the 
double precision is used, we get comparable results with the results 
shown in |18j . 



Table 9: Example 9. 



5 


MAE 




Iter. 


CPU time (seconds) 


Q"ng 


1.00 x 10" 2 


2.07 x 10" 1 


30 


3 


6.25 x 10" 2 


2.00 x 10" 6 


1.00 x 10" 4 


7.14 x 10 -2 


32 


4 


3.44 x 10" 1 


2.00 x 10" 9 


1.00 x 10" 6 


2.56 x 10~ 2 


54 


5 


3.75 x 1Q- 1 


2.00 x 10" 12 



Example 10. (see [5\l) 

f ( t ) = ! *> - 1 < 10 ' 

10 | 0, elsewhere, 
„ / 50, p = 0, 

F 10 {p) = < i-e-iop iQe-iQp n 

I — ^ p — , V > u. 



Table 10: Example 10. 



5 


MAE 




Iter. 


CPU time (seconds) 




1.00 x 10" 2 


2.09 x 10" 1 


30 


3 


3.13 x 10" 2 


2.00 x 10" 6 


1.00 x 10~ 4 


1.35 x 10~ 2 


32 


4 


9.38 x 10~ 2 


2.00 x 10~ 9 


1.00 x 10" 6 


3.00 x 10~ 3 


54 


4 


2.66 x 10" 1 


2.00 x 10~ 9 
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Figure 10: Example 10: the stability of the approximate solution 



Table 10 shows the stability of the solution obtained by our method 
with respect to the noise level 5. We get an excellent agreement be- 
tween the exact solution and the approximate solution for all the noise 
levels 5 as shown in Figure 10. 



Example 11. (see /3]/, JTffi) 

fn(t) 
Fnip) 



sin(t), < t < 10, 
0, elsewhere, 

1 - e- 10p (>sin(10) + cos(10)) 
1 +p 2 




Figure 11: Example 11: the stability of the approximate solution 

Here the function fn(t) represents the class of periodic functions. It is 
mentioned in [16] that oscillating function can be found with accept- 
able accuracy only for relatively small values of t. In this example the 
best approximation is obtained when the noise level 5 = 10 -6 which is 
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comparable with the results given in [5] and [16]. The reconstruction 
of the function fu(t) for various levels of the noise 5 are given in Fig- 
ure 11. The stability of the proposed method with respect to the noise 
S is shown in Table 11. In this example the parameter m ng is bounded 
by the constant 54 when the noise level S = 10~ 6 and k = 0.3. 



Table 11: Example 11. 



5 


MAE 


m ns 


Iter. 


CPU time (seconds) 




1.00 x 10" 2 


2.47 x 10" 1 


30 


3 


9.38 x 10~ 2 


2.00 x 10" 6 


1.00 x 10~ 4 


4.91 x 10~ 2 


32 


4 


2.50 x 10" 1 


2.00 x 10" 9 


1.00 x 10" 6 


2.46 x 10~ 2 


54 


5 


4.38 x 10" 1 


2.00 x 10" 12 



Example 12. (see 13$, 13$) 

/12(f) = 



tcos(t), 0<t<10, 
0, elsewhere, 

(p 2 - 1) - e~ w P(-l+p 2 + lOp + Wp 3 ) cos(lO) 

(1+p 2 ) 2 
e -wp(2p + 10 + 10p 2 ) sin(10) 
+ (1+p 2 ) 2 ' 




2 4 6 8 10 



Figure 12: Example 12: the stability of the approximate solution 

Here we take an increasing function which oscillates as the variable t 
increases over the interval [0, 10). A poor approximation is obtained 
when the noise level 5 = 10" 2 . Fi gure 12 shows that the exact solution 
can be approximated very well when the noise level 5 = 10 -6 . The 
results of our method are comparable with these of the methods given 
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in [3] and [5]. The stability of our method with respect to the noise 
level is shown in Table 12. 



Table 12: Example 12. 



5 


MAE 




Iter. 


CPU time (seconds) 




1.00 x 10~ 2 


1.37 x 10° 


96 


3 


9.38 x 10~ 2 


2.00 x 10~ 6 


1.00 x 10~ 4 


5.98 x 10- 1 


100 


4 


2.66 x lO" 1 


2.00 x 10~ 9 


1.00 x 10" 6 


2.24 x 1Q- 1 


300 


5 


3.44 x 1Q- 1 


2.00 x 10" 12 



• Example 13. 

f n (t) = e- t , F l3 (p) = -^—. 

l+p 

Here the support of / 13(f) is not compact. From the Laplace transform 
formula one gets 

roc rb roo 

F 13 (p) = / e- t e- pt dt= / e~ ( - 1+p ^dt + / e - {1+p)t dt 

JO JO Jb 

= / fi 3 (t)e~ pt dt + —— := h + h, 
Jo 1+P 

where 5(b) := e~ b . Therefore, I2 can be considered as noise of the data 
^13 0), i-e-, 

Ff 3 (p):=F 13 (p)-5(b), (128) 

where 5(b) := e~ b . In this example the following parameters are used: 
d = 2, k = 10" 1 for 5 = e" 5 and k = 10" 5 for 6 = 10" 8 , 10" 20 
and 10" 30 . Table 13 shows that the error decreases as the parameter b 
increases. The approximate solution obtained by the proposed method 
converges to the function /13(f) as b increases (see Figure 13). 
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Table 13: Example 13. 



b 


MAE 




m$ 


Iter 


CPU time (seconds) 


5 


1.487 x 10" 


-2 


2 


4 


3.125 x 10" 2 


8 


2.183 x 10" 


-4 


2 


4 


3.125 x 10~ 2 


20 


4.517 x 10" 


-9 


2 


4 


3.125 x 10~ 2 


30 


1.205 x 10" 


-13 


2 


4 


3.125 x 10" 2 




Figure 13: Example 13: the stability of the approximate solution 



4 Conclusion 



We have tested the proposed algorithm on the wide class of examples con- 
sidered in the literature. Using the rule (i63j) and the stopping rule (j87|) . 
the number of terms in representation (fT3|) . the discrete data F$(pj), j = 
0, 1,2, . . . ,m, and regularization parameter a ng , which are used in comput- 
ing the approximation f^(t) (see {73])) of the unknown function f(t), are 
obtained automatically. Our numerical experiments show that the computa- 
tion time (CPU time) for approximating the function f(t) is small, namely 
CPU time < 1 seconds, and the proposed iterative scheme and the proposed 
adaptive stopping rule yield stable solution with respect to the noise level 5. 
The proposed method also works for / without compact support as shown 
in Example 13. Moreover, in the proposed method we only use a simple 
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representation (|73p which is based on the kernel of the Laplace transform 
integral, so it can be easily implemented numerically. 
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